q— title: “Experimental Data Analysis” author: “Savannah Weaver” date: “July 2021” output: pdf_document: toc: TRUE —

Packages

Data

This data was collected in the Spring of 2021 in conjuction with a study carried out in Cal Poly’s Herpetology class. Some lizards measured for that primary study were kept to observe physiological changes in response to different climate treatments.

Morphometrics & Hydration

Treatment Groups

variables: - individual lizard ID - temp_tmt_C = temperature treatment - humidity_tmt_percent = humidity treatment (high/low, not actually %) - trial_number = which set of lizards that individual was from - conclusion = how that individual’s experiment ended (died, canceled, or complete)

tmts <- read.csv("./data/exp_tmt_assignment.csv")

Capture Data

variables: - date = date of capture & baseline measurements - individual lizard ID - mass_g = mass in grams - hematocrit_percent = % of blood sample that’s red blood cells - osmolality_mmol_kg = concentration of solutes in blood plasma - type = when the measurements were taken along the course of the experiment (all on capture day)

capture_hydration <- read.csv("./exported_data/capture_hydration.csv",
                             na.strings=c("","NA") # fix empty cells
                             ) %>%
  mutate(# correctly format date-only variable
         date = as.Date(date, format = "%Y-%m-%d")
         ) %>%
  # select only relevant variables
  dplyr::select(date, individual_ID, 
                mass_g, hematocrit_percent, osmolality_mmol_kg
                ) %>%
  dplyr::filter(individual_ID %in% tmts$individual_ID) %>%
  mutate(type = as.factor("capture"))
summary(capture_hydration)
##       date            individual_ID        mass_g      hematocrit_percent
##  Min.   :2021-04-19   Min.   : 31.00   Min.   : 8.20   Min.   :16.00     
##  1st Qu.:2021-04-26   1st Qu.: 57.00   1st Qu.:11.10   1st Qu.:32.75     
##  Median :2021-04-26   Median : 78.00   Median :12.65   Median :36.00     
##  Mean   :2021-04-29   Mean   : 77.46   Mean   :12.18   Mean   :36.08     
##  3rd Qu.:2021-05-03   3rd Qu.: 98.25   3rd Qu.:13.32   3rd Qu.:39.00     
##  Max.   :2021-05-10   Max.   :122.00   Max.   :15.00   Max.   :54.00     
##  osmolality_mmol_kg      type   
##  Min.   :319.0      capture:52  
##  1st Qu.:354.2                  
##  Median :373.0                  
##  Mean   :373.7                  
##  3rd Qu.:392.2                  
##  Max.   :423.0

extract SVL data separately from capture data:

SVL <- read.csv("./exported_data/capture_hydration.csv",
                             na.strings=c("","NA") # fix empty cells
                             ) %>%
  dplyr::select(individual_ID, SVL_mm) %>%
  dplyr::filter(individual_ID %in% tmts$individual_ID)
summary(SVL)
##  individual_ID        SVL_mm     
##  Min.   : 31.00   Min.   :59.00  
##  1st Qu.: 57.00   1st Qu.:65.00  
##  Median : 78.00   Median :68.00  
##  Mean   : 77.46   Mean   :67.62  
##  3rd Qu.: 98.25   3rd Qu.:70.00  
##  Max.   :122.00   Max.   :73.00

extract capture CEWL cloacal temperature separately:

cap_CT <- read.csv("./exported_data/capture_hydration.csv",
                             na.strings=c("","NA") # fix empty cells
                             ) %>%
  dplyr::select(individual_ID, cloacal_temp_C) %>%
  dplyr::filter(individual_ID %in% tmts$individual_ID)
summary(cap_CT)
##  individual_ID    cloacal_temp_C 
##  Min.   : 31.00   Min.   :20.00  
##  1st Qu.: 57.00   1st Qu.:22.00  
##  Median : 78.00   Median :24.00  
##  Mean   : 77.46   Mean   :23.68  
##  3rd Qu.: 98.25   3rd Qu.:25.00  
##  Max.   :122.00   Max.   :28.00  
##                   NA's   :2

Experiment Data

variables: - date = date of measurements - individual lizard ID - mass_g = mass in grams - hematocrit_percent = % of blood sample that’s red blood cells - osmolality_mmol_kg = concentration of solutes in blood plasma (mean of 1-3 replicates) - type = when the measurements were taken along the course of the experiment (either during experimental treatment or after rehab)

exp_dat <- read.csv("./data/experimental_data.csv",
                             na.strings=c("","NA") # fix empty cells
                             ) %>%
                # format date
  dplyr::mutate(date = as.Date(date, format = "%m/%d/%y"),
                type = as.factor(type)
                ) %>%
  # select only variables to be analyzed
  dplyr::select(date, individual_ID, mass_g, 
                hematocrit_percent, type, 
                osmolality_mmol_kg = osmolality_mmol_kg_replicate_mean)
summary(exp_dat)
##       date            individual_ID        mass_g       hematocrit_percent
##  Min.   :2021-04-21   Min.   : 31.00   Min.   : 6.700   Min.   :12.0      
##  1st Qu.:2021-04-28   1st Qu.: 51.25   1st Qu.: 9.875   1st Qu.:23.0      
##  Median :2021-05-07   Median : 87.50   Median :11.250   Median :28.0      
##  Mean   :2021-05-06   Mean   : 77.85   Mean   :11.076   Mean   :27.8      
##  3rd Qu.:2021-05-14   3rd Qu.:101.25   3rd Qu.:12.225   3rd Qu.:33.0      
##  Max.   :2021-05-20   Max.   :122.00   Max.   :14.700   Max.   :43.0      
##                                                         NA's   :19        
##     type    osmolality_mmol_kg
##  exp  :98   Min.   :298.0     
##  rehab:34   1st Qu.:342.0     
##             Median :355.0     
##             Mean   :360.1     
##             3rd Qu.:374.8     
##             Max.   :441.0     
##             NA's   :22

Join Dataframes

Now, attach all the dataframes, only use individuals whose treatment was completed, and add a “day” variable for what day of treatment each lizard/observation was on. I also calculate SMI using the equation created in capture_analysis.

all_dat <- exp_dat %>%
  # join data
  rbind(capture_hydration) %>%
  # add tmt group info
  left_join(tmts, by = "individual_ID") %>%
  dplyr::select(-notes) %>%
  # add SVL value for each obs of each indiv.
  # for computing BCI and scaled mass indices
  left_join(SVL, by = "individual_ID") %>%
  # only use completed experiment runs
  dplyr::filter(conclusion == "complete") %>%
  group_by(individual_ID) %>%
  # reformat a lot of variables
  mutate(capture_date = min(date),
         day = as.numeric(date - capture_date),
         humidity_tmt_percent = as.factor(humidity_tmt_percent),
         individual_ID = as.factor(individual_ID),
         temp_tmt_C = as.factor(temp_tmt_C),
         trial_number = as.factor(trial_number),
         conclusion = as.factor(conclusion),
         SMI = mass_g * ((65.02158/SVL_mm) ^ (3.09059/sqrt(0.8944)))
         )

summary(all_dat)
##       date            individual_ID     mass_g      hematocrit_percent
##  Min.   :2021-04-19   37     :  6   Min.   : 6.70   Min.   :12.00     
##  1st Qu.:2021-04-30   39     :  6   1st Qu.:10.20   1st Qu.:24.00     
##  Median :2021-05-07   40     :  6   Median :11.50   Median :30.00     
##  Mean   :2021-05-06   49     :  6   Mean   :11.27   Mean   :29.58     
##  3rd Qu.:2021-05-13   52     :  6   3rd Qu.:12.60   3rd Qu.:35.00     
##  Max.   :2021-05-20   47     :  5   Max.   :15.00   Max.   :54.00     
##                       (Other):116                   NA's   :12        
##       type    osmolality_mmol_kg temp_tmt_C humidity_tmt_percent trial_number
##  exp    :82   Min.   :298.0      25:151     dry  :74             1:35        
##  rehab  :34   1st Qu.:342.8                 humid:77             2:24        
##  capture:35   Median :359.0                                      3:44        
##               Mean   :362.5                                      4:48        
##               3rd Qu.:379.0                                                  
##               Max.   :441.0                                                  
##               NA's   :15                                                     
##     conclusion      SVL_mm       capture_date             day        
##  complete:151   Min.   :59.00   Min.   :2021-04-19   Min.   : 0.000  
##                 1st Qu.:66.00   1st Qu.:2021-04-26   1st Qu.: 2.000  
##                 Median :68.00   Median :2021-05-03   Median : 4.000  
##                 Mean   :67.45   Mean   :2021-04-30   Mean   : 5.424  
##                 3rd Qu.:70.00   3rd Qu.:2021-05-10   3rd Qu.: 9.000  
##                 Max.   :73.00   Max.   :2021-05-10   Max.   :11.000  
##                                                                      
##       SMI        
##  Min.   : 7.343  
##  1st Qu.: 8.990  
##  Median :10.011  
##  Mean   : 9.983  
##  3rd Qu.:10.751  
##  Max.   :13.970  
## 
unique(all_dat$individual_ID)
##  [1] 40  52  37  39  47  49  80  66  54  61  74  73  92  91  95  88  93  96  98 
## [20] 89  99  81  97  104 108 122 118 109 113 105 114 101 117 102 103
## 35 Levels: 37 39 40 47 49 52 54 61 66 73 74 80 81 88 89 91 92 93 95 96 ... 122

re-order some factors:

all_dat$humidity_tmt_percent <- factor(all_dat$humidity_tmt_percent,
                                       levels = c("humid", "dry"),
                                       labels = c("Humid", "Dry"))

make a sub-dataframe without rehab data to prevent any mix-ups:

all_dat_no_rehab <- all_dat %>%
  dplyr::filter(type != "rehab")

Checks

Dates:

# check that capture dates are valid
unique(all_dat$capture_date)
## [1] "2021-04-19" "2021-04-26" "2021-05-03" "2021-05-10"

Check that each lizard only has an accurate number of measurements.

all_dat %>%
  group_by(individual_ID, type) %>%
  summarise(n = n()) %>%
  arrange(type)
## `summarise()` regrouping output by 'individual_ID' (override with `.groups` argument)
## # A tibble: 104 x 3
## # Groups:   individual_ID [35]
##    individual_ID type      n
##    <fct>         <fct> <int>
##  1 37            exp       4
##  2 39            exp       4
##  3 40            exp       4
##  4 47            exp       4
##  5 49            exp       4
##  6 52            exp       4
##  7 54            exp       2
##  8 61            exp       2
##  9 66            exp       2
## 10 73            exp       2
## # … with 94 more rows

That all looks good, experimental measurements are either 4 (first trial) or 2 (other trials). I am excluding lizards that died in treatment from the analysis.

CEWL

Capture CEWL

variables: - date = date of capture & baseline measurements - individual lizard ID - region = which body area the measurement was taken from - TEWL_g_m2h = evaporative water loss - cloacal_temp_C = taken at measurement; influences CEWL

cap_CEWL <- read.csv("./exported_data/capture_CEWL.csv") %>%
  dplyr::select(date, individual_ID, region, TEWL_g_m2h) %>%
  mutate(#individual_ID = as.factor(individual_ID), # do later
         date = as.Date(date, format = "%Y-%m-%d"),
         region = as.factor(region),
         day = as.factor("before"),
         n_day = 0
         ) %>%
  dplyr::filter(individual_ID %in% all_dat$individual_ID) %>%
  left_join(cap_CT, by = 'individual_ID')
summary(cap_CEWL)
##       date            individual_ID     region     TEWL_g_m2h        day     
##  Min.   :2021-04-19   Min.   : 37.00   dewl:32   Min.   : 7.48   before:163  
##  1st Qu.:2021-04-26   1st Qu.: 73.00   dors:33   1st Qu.:20.54               
##  Median :2021-05-03   Median : 95.00   head:33   Median :27.43               
##  Mean   :2021-05-02   Mean   : 87.46   mite:32   Mean   :29.30               
##  3rd Qu.:2021-05-10   3rd Qu.:104.00   vent:33   3rd Qu.:36.91               
##  Max.   :2021-05-10   Max.   :122.00             Max.   :62.94               
##      n_day   cloacal_temp_C 
##  Min.   :0   Min.   :20.00  
##  1st Qu.:0   1st Qu.:22.00  
##  Median :0   Median :24.00  
##  Mean   :0   Mean   :23.84  
##  3rd Qu.:0   3rd Qu.:25.00  
##  Max.   :0   Max.   :28.00

Post-Experiment CEWL

In the future, I could automate this like I did for the HOBO data.

Load in each of the post-rehab datafiles:

# trial 1
CEWL_t1 <- read.csv("./data/post_exp_CEWL/4-28-21-CEWL.csv", # filename
                          na.strings=c("","NA")) %>% # fix empty cells
  # rename and select the pertinent variables/cols
  # I have to do this for each one 
  # so they all have the same number of columns for joining
  dplyr::select(date = Date, 
                Status,
                ID = Comments,
                TEWL_g_m2h = TEWL..g..m2h.. # rename
                )

# trial 2
CEWL_t2 <- read.csv("./data/post_exp_CEWL/5-4-21-CEWL.csv",
                          na.strings=c("","NA")) %>%
  dplyr::select(date = Date, 
                Status,
                ID = Comments,
                TEWL_g_m2h = TEWL..g..m2h..
                )

# trial 3
CEWL_t3 <- read.csv("./data/post_exp_CEWL/5-11-21-CEWL.csv",
                          na.strings=c("","NA")) %>%
  dplyr::select(date = Date, 
                Status,
                ID = Comments,
                TEWL_g_m2h = TEWL..g..m2h..
                )

# trial 4
CEWL_t4 <- read.csv("./data/post_exp_CEWL/5-18-21-CEWL.csv",
                          na.strings=c("","NA")) %>%
  dplyr::select(date = Date, 
                Status,
                ID = Comments,
                TEWL_g_m2h = TEWL..g..m2h..
                )

Load in cloacal temperatures:

exp_CT <- read.csv("./data/post_exp_CEWL_cloacal_temps.csv") %>%
  mutate(date = as.Date(date, format = "%Y/%m/%d")) %>%
  dplyr::select(-time)
summary(exp_CT)
##       date            individual_ID    cloacal_temp_C
##  Min.   :2021-04-28   Min.   : 37.00   Min.   :19.0  
##  1st Qu.:2021-05-04   1st Qu.: 69.50   1st Qu.:21.0  
##  Median :2021-05-11   Median : 93.00   Median :23.0  
##  Mean   :2021-05-09   Mean   : 85.91   Mean   :22.4  
##  3rd Qu.:2021-05-18   3rd Qu.:103.50   3rd Qu.:23.0  
##  Max.   :2021-05-18   Max.   :122.00   Max.   :26.0

Join Dataframes

define not-in function:

`%nin%` = Negate(`%in%`)

Merge all post-experiment CEWL, add cloacal temperature, add capture CEWL:

# merge all CEWL datafiles & reformat
CEWL <- CEWL_t1 %>% # trial 1
  rbind(., CEWL_t2, # trial 2
        CEWL_t3, # trial 3
        CEWL_t4 # trial 4
        ) %>%
  # remove any unsuccessful measurements
  dplyr::filter(Status == "Normal") %>%
  # extract individual_ID and region separately from the "ID" variable
  separate(ID, c("individual_ID", "region")) %>%
  # reformat data
  dplyr::mutate(# reformat date
                date = as.Date(date, format = "%m/%d/%y"),
                # format individual ID
                individual_ID = as.integer(individual_ID),
                # set body region as a factor variable after getting only the consistent characters due to typos
                region = as.factor(substring(region, 1, 4)),
                # add when measurement taken
                day = as.factor("after"),
                n_day = 1 # technically day 8/9, just to help with figures
                ) %>%
  # remove cols not relevant to stats
  dplyr::select(-Status) %>%
  # remove any rows with missing values
  # none actually needed to be removed
  dplyr::filter(complete.cases(.)) %>%
  # add cloacal temperatures
  left_join(exp_CT, by = c("date", "individual_ID")) %>%
  # now matching dataframes, add capture CEWL data
  rbind(cap_CEWL) %>%
  # add tmt assignments
  left_join(tmts, by = "individual_ID") %>%
  mutate(humidity_tmt_percent = as.factor(humidity_tmt_percent),
         individual_ID = as.factor(individual_ID),
         conclusion = as.factor(conclusion),
         trial_number = as.factor(trial_number)
         ) %>%
  # lizards 49 & 80 are missing pre-exp CEWL, so remove them
  dplyr::filter((individual_ID %nin% c('49', '80')))
# every lizard should have 10 measurements
summary(CEWL)
##       date            individual_ID  region     TEWL_g_m2h         day     
##  Min.   :2021-04-19   37     : 10   dewl:65   Min.   :  4.60   after :163  
##  1st Qu.:2021-05-03   39     : 10   dors:65   1st Qu.: 20.09   before:163  
##  Median :2021-05-10   40     : 10   head:66   Median : 27.18               
##  Mean   :2021-05-06   47     : 10   mite:64   Mean   : 30.69               
##  3rd Qu.:2021-05-11   52     : 10   vent:66   3rd Qu.: 38.72               
##  Max.   :2021-05-18   54     : 10             Max.   :106.38               
##                       (Other):266                                          
##      n_day     cloacal_temp_C    temp_tmt_C humidity_tmt_percent trial_number
##  Min.   :0.0   Min.   :19.00   Min.   :25   dry  :158            1: 50       
##  1st Qu.:0.0   1st Qu.:21.00   1st Qu.:25   humid:168            2: 48       
##  Median :0.5   Median :23.00   Median :25                        3:110       
##  Mean   :0.5   Mean   :23.11   Mean   :25                        4:118       
##  3rd Qu.:1.0   3rd Qu.:24.75   3rd Qu.:25                                    
##  Max.   :1.0   Max.   :28.00   Max.   :25                                    
##                                                                              
##     conclusion     notes          
##  complete:326   Length:326        
##                 Class :character  
##                 Mode  :character  
##                                   
##                                   
##                                   
## 

Check that data looks correct:

CEWL %>%
  group_by(individual_ID, day) %>%
  summarise(n = n()) %>%
  arrange(individual_ID, n)
## `summarise()` regrouping output by 'individual_ID' (override with `.groups` argument)
## # A tibble: 66 x 3
## # Groups:   individual_ID [33]
##    individual_ID day        n
##    <fct>         <fct>  <int>
##  1 37            after      5
##  2 37            before     5
##  3 39            after      5
##  4 39            before     5
##  5 40            after      5
##  6 40            before     5
##  7 47            after      5
##  8 47            before     5
##  9 52            after      5
## 10 52            before     5
## # … with 56 more rows

Everything looks great! (after removing the observations for the two lizards with missing pre-experiment CEWL measurements.)

Before/after aren’t perfectly even because sometimes we were unable to get the AquaFlux to equilibrate and take a measurement.

Finally, make a small edit so the regions are spelled out completely. This requires reordering factor levels:

CEWL$region <- factor(CEWL$region, 
                      levels = c("dors", "vent", "head", "dewl", "mite"),
                      labels = c("Dorsum", "Ventrum", "Head", 
                                 "Dewlap", "Mite Patch")
                  )
CEWL$humidity_tmt_percent <- factor(CEWL$humidity_tmt_percent,
                                       levels = c("humid", "dry"),
                                       labels = c("Humid", "Dry"))
CEWL$day <- factor(CEWL$day,
                                       levels = c("before", "after"),
                                       labels = c("Before", "After"))
summary(CEWL)
##       date            individual_ID        region     TEWL_g_m2h    
##  Min.   :2021-04-19   37     : 10   Dorsum    :65   Min.   :  4.60  
##  1st Qu.:2021-05-03   39     : 10   Ventrum   :66   1st Qu.: 20.09  
##  Median :2021-05-10   40     : 10   Head      :66   Median : 27.18  
##  Mean   :2021-05-06   47     : 10   Dewlap    :65   Mean   : 30.69  
##  3rd Qu.:2021-05-11   52     : 10   Mite Patch:64   3rd Qu.: 38.72  
##  Max.   :2021-05-18   54     : 10                   Max.   :106.38  
##                       (Other):266                                   
##      day          n_day     cloacal_temp_C    temp_tmt_C humidity_tmt_percent
##  Before:163   Min.   :0.0   Min.   :19.00   Min.   :25   Humid:168           
##  After :163   1st Qu.:0.0   1st Qu.:21.00   1st Qu.:25   Dry  :158           
##               Median :0.5   Median :23.00   Median :25                       
##               Mean   :0.5   Mean   :23.11   Mean   :25                       
##               3rd Qu.:1.0   3rd Qu.:24.75   3rd Qu.:25                       
##               Max.   :1.0   Max.   :28.00   Max.   :25                       
##                                                                              
##  trial_number    conclusion     notes          
##  1: 50        complete:326   Length:326        
##  2: 48                       Class :character  
##  3:110                       Mode  :character  
##  4:118                                         
##                                                
##                                                
## 

Export Data Frames for Power Analyses

#write.csv(all_dat, "exported_data/exp_effects_hydration.csv")
#write.csv(CEWL, "exported_data/exp_effects_CEWL.csv")

Data Distributions

Histograms

Mass

all_dat %>%
  ggplot(., aes(x = mass_g)) +
  geom_histogram(color = "black", fill="steelblue", bins=50) + 
  theme_classic() +
  xlab("Mass (g)") + 
  ylab("Count")

simple.eda(all_dat$mass_g)

shapiro.test(all_dat$mass_g)
## 
##  Shapiro-Wilk normality test
## 
## data:  all_dat$mass_g
## W = 0.97747, p-value = 0.01396

Mass distribution not normal, skewed to the left.

Scaled Mass Index

all_dat %>%
  ggplot(., aes(x = SMI)) +
  geom_histogram(color = "black", fill="steelblue", bins=50) + 
  theme_classic() +
  xlab("Scaled Mass Index (g)") + 
  ylab("Count")

simple.eda(all_dat$SMI)

shapiro.test(all_dat$SMI)
## 
##  Shapiro-Wilk normality test
## 
## data:  all_dat$SMI
## W = 0.99012, p-value = 0.3712

Hematocrit

all_dat %>%
  ggplot(., aes(x = hematocrit_percent)) +
  geom_histogram(color = "black", fill="steelblue", bins=50) + 
  theme_classic() +
  xlab("Hematocrit (%)") + 
  ylab("Count")
## Warning: Removed 12 rows containing non-finite values (stat_bin).

simple.eda(all_dat$hematocrit_percent)

shapiro.test(all_dat$hematocrit_percent)
## 
##  Shapiro-Wilk normality test
## 
## data:  all_dat$hematocrit_percent
## W = 0.98984, p-value = 0.4089

Visually, looks slightly skewed to the right, but statistically, the distribution of hematocrit is normal.

Osmolality

all_dat %>%
  ggplot(., aes(x = osmolality_mmol_kg)) +
  geom_histogram(color = "black", fill="steelblue", bins=50) + 
  theme_classic() +
  xlab("Osmolality (mmol / kg)") + 
  ylab("Count")
## Warning: Removed 15 rows containing non-finite values (stat_bin).

simple.eda(all_dat$osmolality_mmol_kg)

shapiro.test(all_dat$osmolality_mmol_kg)
## 
##  Shapiro-Wilk normality test
## 
## data:  all_dat$osmolality_mmol_kg
## W = 0.98331, p-value = 0.09544

Visually, looks slightly skewed to the right, but statistically, the distribution of osmolality is normal.

Data Removal

test the effect of removing the day 2/6 data, for consistency’s sake

trimmed <- all_dat_no_rehab %>%
  dplyr::filter(day %in% c(0, 4, 8, 9))
summary(trimmed)
##       date            individual_ID     mass_g      hematocrit_percent
##  Min.   :2021-04-19   37     : 3    Min.   : 6.70   Min.   :12.00     
##  1st Qu.:2021-04-30   39     : 3    1st Qu.:10.20   1st Qu.:27.00     
##  Median :2021-05-07   40     : 3    Median :11.60   Median :32.00     
##  Mean   :2021-05-05   47     : 3    Mean   :11.44   Mean   :31.39     
##  3rd Qu.:2021-05-11   49     : 3    3rd Qu.:12.80   3rd Qu.:36.00     
##  Max.   :2021-05-18   52     : 3    Max.   :15.00   Max.   :54.00     
##                       (Other):87                                      
##       type    osmolality_mmol_kg temp_tmt_C humidity_tmt_percent trial_number
##  exp    :70   Min.   :313.0      25:105     Humid:54             1:18        
##  rehab  : 0   1st Qu.:341.0                 Dry  :51             2:18        
##  capture:35   Median :356.5                                      3:33        
##               Mean   :361.0                                      4:36        
##               3rd Qu.:376.8                                                  
##               Max.   :426.0                                                  
##               NA's   :3                                                      
##     conclusion      SVL_mm      capture_date             day       
##  complete:105   Min.   :59.0   Min.   :2021-04-19   Min.   :0.000  
##                 1st Qu.:66.0   1st Qu.:2021-04-26   1st Qu.:0.000  
##                 Median :68.0   Median :2021-05-03   Median :4.000  
##                 Mean   :67.6   Mean   :2021-05-01   Mean   :4.057  
##                 3rd Qu.:70.0   3rd Qu.:2021-05-10   3rd Qu.:8.000  
##                 Max.   :73.0   Max.   :2021-05-10   Max.   :9.000  
##                                                                    
##       SMI        
##  Min.   : 7.604  
##  1st Qu.: 9.142  
##  Median :10.052  
##  Mean   :10.061  
##  3rd Qu.:10.884  
##  Max.   :13.970  
## 

re-order: NOTE: same difference to make categorical or continuous once day 2/6 are removed

trimmed$day <- factor(trimmed$day,
                      levels = c("0", "4", "8", "9"),
                      labels = c("Before", "Mid", "After", "After"))

Basic Figures & Models

Figure Overlays

all_dat_trimmed_means <- trimmed %>%
  group_by(humidity_tmt_percent, day) %>%
  summarise(SMI_mean = mean(SMI),
            hct_mean = mean(hematocrit_percent))
## `summarise()` regrouping output by 'humidity_tmt_percent' (override with `.groups` argument)
all_dat_trimmed_mean_osml <- trimmed %>%
  dplyr::filter(complete.cases(osmolality_mmol_kg)) %>%
  group_by(humidity_tmt_percent, day) %>%
  summarise(osml_mean = mean(osmolality_mmol_kg))
## `summarise()` regrouping output by 'humidity_tmt_percent' (override with `.groups` argument)

Mass ~ Time

I won’t be using this, SMI is more applicable.

Just look at plot:

all_dat_no_rehab %>% 
  ggplot(data = .) + 
  geom_point(aes(x = day,
                 y = mass_g, 
                 color = humidity_tmt_percent
                 ), 
             size = 1, 
             alpha = 0.6) +
  stat_smooth(aes(x = day, 
                  y = mass_g, 
                  color = humidity_tmt_percent
                  ), 
              formula = y ~ x, 
              method = "lm", 
              se = F, 
              size = 1.6, 
              alpha = 1 ) + 
  theme_classic() + 
  scale_x_continuous(breaks = c(0, 2, 4, 6, 8)) +
  
  scale_color_brewer(palette = "Set2",
                     name = "Humidity Treatment") +
  xlab("Day") + 
  ylab("Mass (g)") + 
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 10),
        legend.text.align = 0
)

SMI ~ Time

plot over course of experiment:

ggplot() + 
  geom_point(data = trimmed,
             aes(x = day,
                 y = SMI, 
                 color = humidity_tmt_percent
                 ), 
             size = 1, 
             alpha = 0.6) +
  geom_line(data = trimmed,
            aes(x = day,
                y = SMI,
                group = individual_ID,
                color = humidity_tmt_percent),
            alpha = 0.2) +
  geom_line(data = all_dat_trimmed_means,
            aes(x = day,
                y = SMI_mean,
                group = humidity_tmt_percent,
                color = humidity_tmt_percent),
            size = 1.6,
            alpha = 1) +
  theme_classic() + 
  scale_color_brewer(palette = "Set2",
                     name = "Humidity Treatment") +
  xlab("") + 
  ylab("Body Condition (g)") + 
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 18),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 18),
        legend.text.align = 0,
        legend.position = "none"
) -> tmt_effects_SMI
tmt_effects_SMI

# export figure
ggsave(filename = "tmt_effects_SMI.jpeg",
       plot = tmt_effects_SMI,
       path = "./final_figures",
       device = "jpeg",
       dpi = 1200,
       width = 5, height = 4)

Check whether means started out different:

SMI_diff_lm <- all_dat_no_rehab %>%
  dplyr::filter(day == 0) %>%
  lm(data = ., SMI ~ humidity_tmt_percent)
summary(SMI_diff_lm)
## 
## Call:
## lm(formula = SMI ~ humidity_tmt_percent, data = .)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9405 -0.7429 -0.0401  0.7385  3.2183 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              10.7517     0.2811   38.25   <2e-16 ***
## humidity_tmt_percentDry  -0.2904     0.4033   -0.72    0.476    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.192 on 33 degrees of freedom
## Multiple R-squared:  0.01547,    Adjusted R-squared:  -0.01436 
## F-statistic: 0.5187 on 1 and 33 DF,  p-value: 0.4765

NOT significantly different, which is good.

Check whether means ended differently:

SMI_diff_lm_end <- all_dat_no_rehab %>%
  dplyr::filter(day %in% c(8, 9)) %>%
  lm(data = ., SMI ~ humidity_tmt_percent)
summary(SMI_diff_lm_end)
## 
## Call:
## lm(formula = SMI ~ humidity_tmt_percent, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.85144 -0.69873 -0.07453  0.80895  2.31088 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               9.8852     0.2585  38.246   <2e-16 ***
## humidity_tmt_percentDry  -0.6287     0.3709  -1.695   0.0994 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.097 on 33 degrees of freedom
## Multiple R-squared:  0.08012,    Adjusted R-squared:  0.05224 
## F-statistic: 2.874 on 1 and 33 DF,  p-value: 0.09943

mixed model:

SMI_mod <- lme4::lmer(data = all_dat_no_rehab,
               SMI ~ day*humidity_tmt_percent +
               (1|trial_number))
summary(SMI_mod)
## Linear mixed model fit by REML ['lmerMod']
## Formula: SMI ~ day * humidity_tmt_percent + (1 | trial_number)
##    Data: all_dat_no_rehab
## 
## REML criterion at convergence: 369.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.84407 -0.68219 -0.04622  0.62984  2.74822 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  trial_number (Intercept) 0.01241  0.1114  
##  Residual                 1.27504  1.1292  
## Number of obs: 117, groups:  trial_number, 4
## 
## Fixed effects:
##                             Estimate Std. Error t value
## (Intercept)                 10.79499    0.24058  44.871
## day                         -0.10598    0.04510  -2.350
## humidity_tmt_percentDry     -0.40660    0.33497  -1.214
## day:humidity_tmt_percentDry -0.04498    0.06463  -0.696
## 
## Correlation of Fixed Effects:
##             (Intr) day    hmd__D
## day         -0.759              
## hmdty_tmt_D -0.678  0.545       
## dy:hmdty__D  0.530 -0.698 -0.782
drop1(SMI_mod) # this was more intuitive using lme4::lmer but I switched to be able to get p-values
## Single term deletions
## 
## Model:
## SMI ~ day * humidity_tmt_percent + (1 | trial_number)
##                          npar    AIC
## <none>                        369.24
## day:humidity_tmt_percent    1 367.74
# drop interaction term
SMI_mod2 <- lmerTest::lmer(data = all_dat_no_rehab,
               SMI ~ day + humidity_tmt_percent +
               (1|trial_number))
summary(SMI_mod2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: SMI ~ day + humidity_tmt_percent + (1 | trial_number)
##    Data: all_dat_no_rehab
## 
## REML criterion at convergence: 366.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.84939 -0.64254 -0.07319  0.69627  2.67483 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  trial_number (Intercept) 0.01262  0.1123  
##  Residual                 1.26911  1.1265  
## Number of obs: 117, groups:  trial_number, 4
## 
## Fixed effects:
##                          Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)              10.88366    0.20374  20.15193  53.420  < 2e-16 ***
## day                      -0.12788    0.03223 110.92424  -3.967 0.000129 ***
## humidity_tmt_percentDry  -0.58888    0.20841 111.10846  -2.826 0.005597 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) day   
## day         -0.641       
## hmdty_tmt_D -0.498  0.000
#drop1(SMI_mod2)
SMI_mod2t <- lmerTest::lmer(data = trimmed,
               SMI ~ day + humidity_tmt_percent +
               (1|trial_number))
## boundary (singular) fit: see ?isSingular
summary(SMI_mod2t)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: SMI ~ day + humidity_tmt_percent + (1 | trial_number)
##    Data: trimmed
## 
## REML criterion at convergence: 327.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.76998 -0.68801 -0.00593  0.73878  2.75012 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  trial_number (Intercept) 0.000    0.000   
##  Residual                 1.303    1.141   
## Number of obs: 105, groups:  trial_number, 4
## 
## Fixed effects:
##                         Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)              10.8313     0.2212 101.0000  48.964  < 2e-16 ***
## dayMid                   -0.6173     0.2728 101.0000  -2.263 0.025803 *  
## dayAfter                 -1.0308     0.2728 101.0000  -3.778 0.000267 ***
## humidity_tmt_percentDry  -0.4543     0.2229 101.0000  -2.039 0.044091 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) dayMid dyAftr
## dayMid      -0.617              
## dayAfter    -0.617  0.500       
## hmdty_tmt_D -0.489  0.000  0.000
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
#drop1(SMI_mod2)

SMI is best predicted by day and treatment, but not including their interaction.

write.csv(broom.mixed::tidy(SMI_mod2t), 
          "./best models/exp_effects_SMI.csv")

Hct ~ Time

ggplot() + 
  geom_point(data = trimmed, 
             aes(x = day,
                 y = hematocrit_percent, 
                 color = humidity_tmt_percent
                 ), 
             size = 1, 
             alpha = 0.6) +
  geom_line(data = trimmed,
            aes(x = day,
                y = hematocrit_percent,
                group = individual_ID,
                color = humidity_tmt_percent),
            alpha = 0.2) +
  geom_line(data = all_dat_trimmed_means,
            aes(x = day,
                y = hct_mean,
                group = humidity_tmt_percent,
                color = humidity_tmt_percent),
            size = 1.6,
            alpha = 1) +
  theme_classic() + 
  scale_color_brewer(palette = "Set2",
                     name = "Humidity Treatment") +
  xlab("Day") + 
  ylab("Hematocrit (%)") + 
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 18),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.text.align = 0,
        legend.position = "none"
        ) -> tmt_effects_hct
tmt_effects_hct

# export figure
ggsave(filename = "tmt_effects_hct.jpeg",
       plot = tmt_effects_hct,
       path = "./final_figures",
       device = "jpeg",
       dpi = 1200,
       width = 5, height = 4)

this model seemed to work well with indiv as a random factor, but still excluded because it’s probably unnecessary

hct_mod <- all_dat_no_rehab %>%
  dplyr::filter(complete.cases(hematocrit_percent)) %>%
  lme4::lmer(data = .,
               hematocrit_percent ~ day + humidity_tmt_percent +
               (1|trial_number))
summary(hct_mod)
## Linear mixed model fit by REML ['lmerMod']
## Formula: hematocrit_percent ~ day + humidity_tmt_percent + (1 | trial_number)
##    Data: .
## 
## REML criterion at convergence: 765.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2154 -0.6104  0.0919  0.6453  2.7070 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  trial_number (Intercept)  8.339   2.888   
##  Residual                 40.121   6.334   
## Number of obs: 117, groups:  trial_number, 4
## 
## Fixed effects:
##                         Estimate Std. Error t value
## (Intercept)              35.2145     1.8208  19.340
## day                      -1.0534     0.1813  -5.811
## humidity_tmt_percentDry  -0.2733     1.1724  -0.233
## 
## Correlation of Fixed Effects:
##             (Intr) day   
## day         -0.403       
## hmdty_tmt_D -0.314  0.000
drop1(hct_mod)
## Single term deletions
## 
## Model:
## hematocrit_percent ~ day + humidity_tmt_percent + (1 | trial_number)
##                      npar    AIC
## <none>                    778.55
## day                     1 806.68
## humidity_tmt_percent    1 776.61
# drop humidity
hct_mod2 <- all_dat_no_rehab %>%
  dplyr::filter(complete.cases(hematocrit_percent)) %>%
  lmerTest::lmer(data = .,
               hematocrit_percent ~ day +
               (1|trial_number))
summary(hct_mod2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: hematocrit_percent ~ day + (1 | trial_number)
##    Data: .
## 
## REML criterion at convergence: 767.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2511 -0.5977  0.1135  0.6655  2.7398 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  trial_number (Intercept)  8.36    2.891   
##  Residual                 39.78    6.307   
## Number of obs: 117, groups:  trial_number, 4
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  35.0810     1.7278   4.5115  20.304 1.31e-05 ***
## day          -1.0534     0.1805 112.0809  -5.836 5.29e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## day -0.423
drop1(hct_mod2)
## Single term deletions using Satterthwaite's method:
## 
## Model:
## hematocrit_percent ~ day + (1 | trial_number)
##     Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## day 1354.8  1354.8     1 112.08  34.057 5.287e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hct_mod2t <- trimmed %>%
  dplyr::filter(complete.cases(hematocrit_percent)) %>%
  lmerTest::lmer(data = .,
               hematocrit_percent ~ day +
               (1|trial_number))
summary(hct_mod2t)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: hematocrit_percent ~ day + (1 | trial_number)
##    Data: .
## 
## REML criterion at convergence: 676.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.14009 -0.58850  0.00871  0.64072  2.74219 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  trial_number (Intercept)  7.253   2.693   
##  Residual                 37.932   6.159   
## Number of obs: 105, groups:  trial_number, 4
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)   35.592      1.712  4.853  20.794 6.24e-06 ***
## dayMid        -5.771      1.472 98.766  -3.920 0.000163 ***
## dayAfter      -8.057      1.472 98.766  -5.473 3.36e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) dayMid
## dayMid   -0.430       
## dayAfter -0.430  0.500

The model AIC is slightly better without the interaction effect, so I removed that. The effect of humidity could ALSO be dropped, so humidity treatment was not an important factor affecting hematocrit, but how many days lizards were in treatment was. Both treatment groups lost hematocrit at approximately the same rate.

write.csv(broom::tidy(hct_mod2t), 
          "./best models/exp_effects_hct.csv")

Osml ~ Time

ggplot() + 
  geom_point(data = trimmed,
             aes(x = day,
                 y = osmolality_mmol_kg, 
                 color = humidity_tmt_percent
                 ), 
             size = 1, 
             alpha = 0.6) +
  geom_line(data = trimmed,
            aes(x = day,
                y = osmolality_mmol_kg,
                group = individual_ID,
                color = humidity_tmt_percent),
            alpha = 0.2) +
  geom_line(data = all_dat_trimmed_mean_osml,
            aes(x = day,
                y = osml_mean,
                group = humidity_tmt_percent,
                color = humidity_tmt_percent),
            size = 1.6,
            alpha = 1) +
  theme_classic() + 
  scale_color_brewer(palette = "Set2",
                     name = "Humidity Treatment") +
  xlab("") + 
  ylab("Osmolality (mmol / kg)") + 
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 18),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.text.align = 0,
        legend.position = "none"
        ) -> tmt_effects_osml
tmt_effects_osml
## Warning: Removed 3 rows containing missing values (geom_point).
## Warning: Removed 3 row(s) containing missing values (geom_path).

# export figure
ggsave(filename = "tmt_effects_osml.jpeg",
       plot = tmt_effects_osml,
       path = "./final_figures",
       device = "jpeg",
       dpi = 1200,
       width = 5, height = 4)
## Warning: Removed 3 rows containing missing values (geom_point).

## Warning: Removed 3 row(s) containing missing values (geom_path).

singular warning - do NOT include individual ID as a random effect

osml_mod <- all_dat_no_rehab %>%
  dplyr::filter(complete.cases(osmolality_mmol_kg)) %>%
  lmerTest::lmer(data = .,
               osmolality_mmol_kg ~ day * humidity_tmt_percent +
               (1|trial_number))
summary(osml_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: osmolality_mmol_kg ~ day * humidity_tmt_percent + (1 | trial_number)
##    Data: .
## 
## REML criterion at convergence: 1018.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1763 -0.6722 -0.1769  0.6925  2.5285 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  trial_number (Intercept) 302.6    17.40   
##  Residual                 469.3    21.66   
## Number of obs: 114, groups:  trial_number, 4
## 
## Fixed effects:
##                             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)                  369.938      9.804   4.188  37.733 1.83e-06 ***
## day                           -1.315      0.877 106.917  -1.499    0.137    
## humidity_tmt_percentDry       -8.152      6.439 106.912  -1.266    0.208    
## day:humidity_tmt_percentDry    1.901      1.266 106.908   1.502    0.136    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) day    hmd__D
## day         -0.357              
## hmdty_tmt_D -0.319  0.542       
## dy:hmdty__D  0.246 -0.691 -0.776
drop1(osml_mod)
## Single term deletions using Satterthwaite's method:
## 
## Model:
## osmolality_mmol_kg ~ day * humidity_tmt_percent + (1 | trial_number)
##                          Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## day:humidity_tmt_percent 1058.1  1058.1     1 106.91  2.2545 0.1362
osml_modt <- trimmed %>%
  dplyr::filter(complete.cases(osmolality_mmol_kg)) %>%
  lmerTest::lmer(data = .,
               osmolality_mmol_kg ~ day * humidity_tmt_percent +
               (1|trial_number))
summary(osml_modt)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: osmolality_mmol_kg ~ day * humidity_tmt_percent + (1 | trial_number)
##    Data: .
## 
## REML criterion at convergence: 872.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.66849 -0.79364  0.00361  0.51469  2.41228 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  trial_number (Intercept) 313.5    17.71   
##  Residual                 393.8    19.84   
## Number of obs: 102, groups:  trial_number, 4
## 
## Fixed effects:
##                                  Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)                       378.509     10.035   4.379  37.719 1.13e-06
## dayMid                            -29.500      6.615  92.951  -4.460 2.29e-05
## dayAfter                          -13.885      6.713  92.956  -2.068   0.0414
## humidity_tmt_percentDry           -10.215      6.713  92.956  -1.522   0.1315
## dayMid:humidity_tmt_percentDry     12.324      9.491  92.951   1.298   0.1974
## dayAfter:humidity_tmt_percentDry   17.783      9.721  92.954   1.829   0.0705
##                                     
## (Intercept)                      ***
## dayMid                           ***
## dayAfter                         *  
## humidity_tmt_percentDry             
## dayMid:humidity_tmt_percentDry      
## dayAfter:humidity_tmt_percentDry .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) dayMid dyAftr hmd__D dM:__D
## dayMid      -0.330                            
## dayAfter    -0.326  0.493                     
## hmdty_tmt_D -0.326  0.493  0.485              
## dyMd:hmd__D  0.230 -0.697 -0.343 -0.707       
## dyAftr:h__D  0.224 -0.340 -0.690 -0.690  0.488
drop1(osml_modt)
## Single term deletions using Satterthwaite's method:
## 
## Model:
## osmolality_mmol_kg ~ day * humidity_tmt_percent + (1 | trial_number)
##                          Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## day:humidity_tmt_percent 1402.8  701.41     2 92.953  1.7812 0.1741

The model seems good as-is.

write.csv(broom::tidy(osml_modt), 
          "./best models/exp_effects_osml.csv")

Change in Osmolality

osml_d0 <- all_dat_no_rehab %>%
  dplyr::filter(day == 0) %>%
  dplyr::select(individual_ID, osml0 = osmolality_mmol_kg, 
                humidity_tmt_percent)
osml_d8 <- all_dat_no_rehab %>%
  dplyr::filter(day %in% c(8,9)) %>%
  dplyr::select(individual_ID, osml89 = osmolality_mmol_kg, 
                humidity_tmt_percent)

osml_diffs <- osml_d0 %>%
  left_join(osml_d8) %>%
  mutate(osml_change = osml89 - osml0)
## Joining, by = c("individual_ID", "humidity_tmt_percent")

boxplot:

osml_diffs %>% 
  ggplot(data = .) + 
  geom_boxplot(aes(x = humidity_tmt_percent, 
                   y = osml_change, 
                   group = humidity_tmt_percent,
                   color = humidity_tmt_percent
                   ), 
               size = 1,
               alpha = 1) + 
  theme_classic() + 
  geom_hline(yintercept = 0) +
  xlab("") + 
  ylab("Osmolality Change (mmol / kg)") + 
  scale_color_brewer(palette = "Set2") +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 18),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 14),
        legend.text.align = 0,
        legend.position = "none"
)
## Warning: Removed 3 rows containing non-finite values (stat_boxplot).

CEWL ~ Before/After

try a boxplot:

CEWL %>% 
  ggplot(data = .) + 
  geom_boxplot(aes(x = humidity_tmt_percent, 
                   y = TEWL_g_m2h, 
                   color = day
                   ), 
               size = 1,
               alpha = 1) + 
  facet_wrap(~region) +
  theme_classic() + 
  xlab("") + 
  ylab("CEWL (g / m^2h)") + 
  scale_color_brewer(palette = "Set2") +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 18),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 14),
        legend.text.align = 0,
        legend.position = "bottom"
)

this is difficult to see changes, I think a line graph would be better…

CEWL %>% 
  ggplot(data = .) + 
  geom_point(aes(x = n_day,
                 y = TEWL_g_m2h, 
                 color = humidity_tmt_percent
                 ), 
             size = 1, 
             alpha = 0.6) +
  geom_line(aes(x = n_day,
                y = TEWL_g_m2h,
                group = individual_ID,
                color = humidity_tmt_percent),
            alpha = 0.2) +
  stat_smooth(aes(x = n_day, 
                  y = TEWL_g_m2h, 
                  color = humidity_tmt_percent
                  ),
              formula = y ~ x, 
              method = "lm", 
              se = F, 
              size = 1.6, 
              alpha = 1) + 
  theme_classic() + 
  scale_color_brewer(palette = "Set2",
                     name = "Humidity Treatment") +
  facet_wrap(~region, ncol = 2) +
  scale_x_continuous(breaks = c(0, 1),
                     labels = c("0" = "before", "1" = "after")
                     ) +
  xlab("") + 
  ylab(bquote('CEWL (g/'*m^2*'h)')) + 
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 10,
                                 angle = 90),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 10),
        legend.text.align = 0,
        legend.position = c(0.75,0.12),
        #legend.justification = c(1, 1)
) -> CEWL_tmt_fig
CEWL_tmt_fig

# export figure
#ggsave(filename = "tmt_effects_CEWL.jpeg",
 #      plot = CEWL_tmt_fig,
  #     path = "./final_figures",
   #    device = "jpeg",
    #   dpi = 1200,
     #  width = 4, height = 6)

I saved the legend separately to make the figure layout better.

CEWL_mod <- CEWL %>%
  lme4::lmer(data = .,
               TEWL_g_m2h ~ day * humidity_tmt_percent * region +
               cloacal_temp_C +
               (1|trial_number/individual_ID))
summary(CEWL_mod)
## Linear mixed model fit by REML ['lmerMod']
## Formula: TEWL_g_m2h ~ day * humidity_tmt_percent * region + cloacal_temp_C +  
##     (1 | trial_number/individual_ID)
##    Data: .
## 
## REML criterion at convergence: 2441.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4361 -0.5707 -0.0822  0.4555  4.1852 
## 
## Random effects:
##  Groups                     Name        Variance Std.Dev.
##  individual_ID:trial_number (Intercept)  38.212   6.182  
##  trial_number               (Intercept)   8.141   2.853  
##  Residual                               123.274  11.103  
## Number of obs: 326, groups:  individual_ID:trial_number, 33; trial_number, 4
## 
## Fixed effects:
##                                                    Estimate Std. Error t value
## (Intercept)                                       -69.91895   11.98454  -5.834
## dayAfter                                           20.55665    3.96057   5.190
## humidity_tmt_percentDry                             0.09390    4.44521   0.021
## regionVentrum                                      14.06471    3.80826   3.693
## regionHead                                          5.20176    3.80826   1.366
## regionDewlap                                        2.76945    3.87231   0.715
## regionMite Patch                                    5.02118    3.80826   1.318
## cloacal_temp_C                                      3.90353    0.48902   7.982
## dayAfter:humidity_tmt_percentDry                  -22.17103    5.52280  -4.014
## dayAfter:regionVentrum                              4.82016    5.43126   0.887
## dayAfter:regionHead                                -7.17102    5.43126  -1.320
## dayAfter:regionDewlap                              -2.93470    5.46212  -0.537
## dayAfter:regionMite Patch                          -6.04716    5.47797  -1.104
## humidity_tmt_percentDry:regionVentrum              -0.34596    5.46920  -0.063
## humidity_tmt_percentDry:regionHead                 -0.03395    5.46920  -0.006
## humidity_tmt_percentDry:regionDewlap                0.90930    5.51399   0.165
## humidity_tmt_percentDry:regionMite Patch            2.90477    5.52066   0.526
## dayAfter:humidity_tmt_percentDry:regionVentrum     -2.37516    7.76641  -0.306
## dayAfter:humidity_tmt_percentDry:regionHead         5.28070    7.76641   0.680
## dayAfter:humidity_tmt_percentDry:regionDewlap       5.66684    7.82338   0.724
## dayAfter:humidity_tmt_percentDry:regionMite Patch   1.47809    7.83657   0.189
## 
## Correlation matrix not shown by default, as p = 21 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
drop1(CEWL_mod)
## Single term deletions
## 
## Model:
## TEWL_g_m2h ~ day * humidity_tmt_percent * region + cloacal_temp_C + 
##     (1 | trial_number/individual_ID)
##                                 npar    AIC
## <none>                               2570.1
## cloacal_temp_C                     1 2628.2
## day:humidity_tmt_percent:region    4 2563.8

Drop triple interaction. I think the day:region standalone would be weird too.

CEWL_mod2 <- CEWL %>%
  lme4::lmer(data = .,
               TEWL_g_m2h ~ 
               humidity_tmt_percent * (day + region) +
               cloacal_temp_C +
               (1|trial_number/individual_ID))
summary(CEWL_mod2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: TEWL_g_m2h ~ humidity_tmt_percent * (day + region) + cloacal_temp_C +  
##     (1 | trial_number/individual_ID)
##    Data: .
## 
## REML criterion at convergence: 2490.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4790 -0.6180 -0.0854  0.4480  4.0055 
## 
## Random effects:
##  Groups                     Name        Variance Std.Dev.
##  individual_ID:trial_number (Intercept)  38.017   6.166  
##  trial_number               (Intercept)   8.677   2.946  
##  Residual                               123.628  11.119  
## Number of obs: 326, groups:  individual_ID:trial_number, 33; trial_number, 4
## 
## Fixed effects:
##                                          Estimate Std. Error t value
## (Intercept)                              -68.6125    11.9022  -5.765
## humidity_tmt_percentDry                   -0.8482     3.7242  -0.228
## dayAfter                                  18.2857     1.9294   9.478
## regionVentrum                             16.5143     2.7189   6.074
## regionHead                                 1.6557     2.7189   0.609
## regionDewlap                               1.3212     2.7199   0.486
## regionMite Patch                           2.0988     2.7414   0.766
## cloacal_temp_C                             3.8933     0.4894   7.956
## humidity_tmt_percentDry:dayAfter         -20.1573     2.4933  -8.084
## humidity_tmt_percentDry:regionVentrum     -1.5730     3.8883  -0.405
## humidity_tmt_percentDry:regionHead         2.5669     3.8883   0.660
## humidity_tmt_percentDry:regionDewlap       3.6718     3.9063   0.940
## humidity_tmt_percentDry:regionMite Patch   3.4656     3.9211   0.884
## 
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
drop1(CEWL_mod2)
## Single term deletions
## 
## Model:
## TEWL_g_m2h ~ humidity_tmt_percent * (day + region) + cloacal_temp_C + 
##     (1 | trial_number/individual_ID)
##                             npar    AIC
## <none>                           2563.5
## cloacal_temp_C                 1 2619.8
## humidity_tmt_percent:day       1 2622.6
## humidity_tmt_percent:region    4 2558.4

We can drop the humidity:region interaction.

CEWL_mod3 <- CEWL %>%
  dplyr::filter(complete.cases(.)) %>%
  lmerTest::lmer(data = .,
               TEWL_g_m2h ~ 
               day*humidity_tmt_percent + region + cloacal_temp_C +
               (1|trial_number/individual_ID))
summary(CEWL_mod3)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: TEWL_g_m2h ~ day * humidity_tmt_percent + region + cloacal_temp_C +  
##     (1 | trial_number/individual_ID)
##    Data: .
## 
## REML criterion at convergence: 2510.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3478 -0.6058 -0.1117  0.4446  3.9319 
## 
## Random effects:
##  Groups                     Name        Variance Std.Dev.
##  individual_ID:trial_number (Intercept)  38.023   6.166  
##  trial_number               (Intercept)   8.872   2.979  
##  Residual                               123.111  11.096  
## Number of obs: 326, groups:  individual_ID:trial_number, 33; trial_number, 4
## 
## Fixed effects:
##                                  Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)                      -69.3758    11.8197 230.9696  -5.869 1.51e-08
## dayAfter                          18.2632     1.9251 304.9541   9.487  < 2e-16
## humidity_tmt_percentDry            0.7605     2.7927  41.1331   0.272   0.7867
## regionVentrum                     15.7651     1.9396 286.0095   8.128 1.33e-14
## regionHead                         2.9138     1.9396 286.0095   1.502   0.1341
## regionDewlap                       3.0977     1.9483 286.0970   1.590   0.1130
## regionMite Patch                   3.7912     1.9565 286.2700   1.938   0.0536
## cloacal_temp_C                     3.8920     0.4886 279.6344   7.966 4.17e-14
## dayAfter:humidity_tmt_percentDry -20.1375     2.4881 288.3429  -8.094 1.63e-14
##                                     
## (Intercept)                      ***
## dayAfter                         ***
## humidity_tmt_percentDry             
## regionVentrum                    ***
## regionHead                          
## regionDewlap                        
## regionMite Patch                 .  
## cloacal_temp_C                   ***
## dayAfter:humidity_tmt_percentDry ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) dyAftr hmd__D rgnVnt regnHd rgnDwl rgnMtP clc__C
## dayAfter    -0.507                                                 
## hmdty_tmt_D  0.005  0.217                                          
## regionVntrm -0.086 -0.004  0.000                                   
## regionHead  -0.086 -0.004  0.000  0.504                            
## regionDewlp -0.098 -0.010 -0.006  0.502  0.502                     
## reginMtPtch -0.108  0.013  0.002  0.500  0.500  0.498              
## clocl_tmp_C -0.972  0.456 -0.122  0.004  0.004  0.017  0.026       
## dyAftr:h__D  0.191 -0.680 -0.418  0.004  0.004  0.017 -0.010 -0.146
drop1(CEWL_mod3)
## Single term deletions using Satterthwaite's method:
## 
## Model:
## TEWL_g_m2h ~ day * humidity_tmt_percent + region + cloacal_temp_C + (1 | trial_number/individual_ID)
##                          Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## region                   9880.0  2470.0     4 286.10  20.063 1.409e-14 ***
## cloacal_temp_C           7812.3  7812.3     1 279.63  63.458 4.165e-14 ***
## day:humidity_tmt_percent 8064.6  8064.6     1 288.34  65.507 1.634e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The model is best with all the parameters currently included in model 3.

write.csv(broom::tidy(CEWL_mod3), 
          "./best models/exp_effects_CEWL.csv")

Change in CEWL

CEWL_before <- CEWL %>%
  dplyr::filter(day == "before") %>%
  dplyr::select(CEWL_before = TEWL_g_m2h,
                humidity_tmt_percent, trial_number,
                individual_ID, region)
CEWL_after <- CEWL %>%
  dplyr::filter(day == "after") %>%
  dplyr::select(CEWL_after = TEWL_g_m2h,
                humidity_tmt_percent, trial_number,
                individual_ID, region)

CEWL_diffs <- CEWL_before %>%
  left_join(CEWL_after, by = c('individual_ID', 'region',
                'humidity_tmt_percent', 'trial_number')) %>%
  mutate(CEWL_diff = CEWL_after - CEWL_before)

plot:

CEWL_diffs %>% 
  ggplot(data = .) + 
  geom_boxplot(aes(x = humidity_tmt_percent, 
                   y = CEWL_diff, 
                   group = humidity_tmt_percent,
                   color = humidity_tmt_percent
                   ), 
               size = 1,
               alpha = 1) + 
  #facet_wrap(~humidity_tmt_percent) +
  theme_classic() + 
  geom_hline(yintercept = 0, lty = 2) +
  xlab("") + 
  ylab("CEWL Change (g/m2h)") + 
  #annotate("text", x = 1.5, y = 45, 
   #        label = "paste(italic(p), \" = 0.0152\")", 
    #       parse = TRUE,
     #      size = 6) +
  #ylim(10, 50) +
  #scale_x_discrete(labels = c("F" = "Female",
   #                           "M" = "Male")) +
  scale_color_brewer(palette = "Set2") +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 18),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 14),
        legend.text.align = 0,
        legend.position = "none"
)

Multi-Figure

ggarrange(tmt_effects_osml, tmt_effects_SMI, tmt_effects_hct,
          ncol = 1, nrow = 3,
          common.legend = TRUE,
          legend = "bottom"
          ) -> tmt_multi_fig
## Warning: Removed 3 rows containing missing values (geom_point).
## Warning: Removed 3 row(s) containing missing values (geom_path).
## Warning: Removed 3 rows containing missing values (geom_point).
## Warning: Removed 3 row(s) containing missing values (geom_path).
tmt_multi_fig

# export figure
ggsave(filename = "tmt_multi_fig.jpeg",
       plot = tmt_multi_fig,
       path = "./final_figures",
       device = "jpeg",
       dpi = 1200,
       width = 6, height = 12)

Rehydration…

Data

First, get only the data for before experiment, after experiment, and after rehab.

summary(all_dat)
##       date            individual_ID     mass_g      hematocrit_percent
##  Min.   :2021-04-19   37     :  6   Min.   : 6.70   Min.   :12.00     
##  1st Qu.:2021-04-30   39     :  6   1st Qu.:10.20   1st Qu.:24.00     
##  Median :2021-05-07   40     :  6   Median :11.50   Median :30.00     
##  Mean   :2021-05-06   49     :  6   Mean   :11.27   Mean   :29.58     
##  3rd Qu.:2021-05-13   52     :  6   3rd Qu.:12.60   3rd Qu.:35.00     
##  Max.   :2021-05-20   47     :  5   Max.   :15.00   Max.   :54.00     
##                       (Other):116                   NA's   :12        
##       type    osmolality_mmol_kg temp_tmt_C humidity_tmt_percent trial_number
##  exp    :82   Min.   :298.0      25:151     Humid:77             1:35        
##  rehab  :34   1st Qu.:342.8                 Dry  :74             2:24        
##  capture:35   Median :359.0                                      3:44        
##               Mean   :362.5                                      4:48        
##               3rd Qu.:379.0                                                  
##               Max.   :441.0                                                  
##               NA's   :15                                                     
##     conclusion      SVL_mm       capture_date             day        
##  complete:151   Min.   :59.00   Min.   :2021-04-19   Min.   : 0.000  
##                 1st Qu.:66.00   1st Qu.:2021-04-26   1st Qu.: 2.000  
##                 Median :68.00   Median :2021-05-03   Median : 4.000  
##                 Mean   :67.45   Mean   :2021-04-30   Mean   : 5.424  
##                 3rd Qu.:70.00   3rd Qu.:2021-05-10   3rd Qu.: 9.000  
##                 Max.   :73.00   Max.   :2021-05-10   Max.   :11.000  
##                                                                      
##       SMI        
##  Min.   : 7.343  
##  1st Qu.: 8.990  
##  Median :10.011  
##  Mean   : 9.983  
##  3rd Qu.:10.751  
##  Max.   :13.970  
## 
rehydrat_dat <- all_dat %>%
  dplyr::filter(day %in% c(0, 8, 9, 10, 11))
rehydrat_dat$day <- factor(rehydrat_dat$day,
                           levels = c(0, 8, 9, 10, 11),
                           labels = c("Before Experiment",
                                      "After Experiment",
                                      "After Experiment",
                                      "After Rehydration",
                                      "After Rehydration"))
summary(rehydrat_dat)
##       date            individual_ID     mass_g      hematocrit_percent
##  Min.   :2021-04-19   37     : 3    Min.   : 6.70   Min.   :14.00     
##  1st Qu.:2021-05-03   39     : 3    1st Qu.:10.20   1st Qu.:23.88     
##  Median :2021-05-10   40     : 3    Median :11.40   Median :29.50     
##  Mean   :2021-05-07   49     : 3    Mean   :11.35   Mean   :29.84     
##  3rd Qu.:2021-05-13   52     : 3    3rd Qu.:12.80   3rd Qu.:35.00     
##  Max.   :2021-05-20   54     : 3    Max.   :15.00   Max.   :54.00     
##                       (Other):86                    NA's   :12        
##       type    osmolality_mmol_kg temp_tmt_C humidity_tmt_percent trial_number
##  exp    :35   Min.   :298.0      25:104     Humid:53             1:17        
##  rehab  :34   1st Qu.:347.0                 Dry  :51             2:18        
##  capture:35   Median :369.0                                      3:33        
##               Mean   :368.5                                      4:36        
##               3rd Qu.:389.0                                                  
##               Max.   :441.0                                                  
##               NA's   :15                                                     
##     conclusion      SVL_mm       capture_date                       day    
##  complete:104   Min.   :59.00   Min.   :2021-04-19   Before Experiment:35  
##                 1st Qu.:66.00   1st Qu.:2021-04-26   After Experiment :35  
##                 Median :68.00   Median :2021-05-03   After Rehydration:34  
##                 Mean   :67.68   Mean   :2021-05-01                         
##                 3rd Qu.:70.00   3rd Qu.:2021-05-10                         
##                 Max.   :73.00   Max.   :2021-05-10                         
##                                                                            
##       SMI        
##  Min.   : 7.343  
##  1st Qu.: 8.965  
##  Median :10.011  
##  Mean   : 9.947  
##  3rd Qu.:10.731  
##  Max.   :13.970  
## 

prepare summary data to overlay on plots:

hydrat_means <- rehydrat_dat %>%
  dplyr::filter(complete.cases(osmolality_mmol_kg, SMI, 
                               hematocrit_percent)) %>%
  group_by(humidity_tmt_percent, day) %>%
  summarise(osml_mean = mean(osmolality_mmol_kg),
            SMI_mean = mean(SMI),
            hct_mean = mean(hematocrit_percent))
## `summarise()` regrouping output by 'humidity_tmt_percent' (override with `.groups` argument)
hydrat_means
## # A tibble: 6 x 5
## # Groups:   humidity_tmt_percent [2]
##   humidity_tmt_percent day               osml_mean SMI_mean hct_mean
##   <fct>                <fct>                 <dbl>    <dbl>    <dbl>
## 1 Humid                Before Experiment      375.    10.8      36.4
## 2 Humid                After Experiment       362.     9.94     27.9
## 3 Humid                After Rehydration      367.     9.72     23.5
## 4 Dry                  Before Experiment      365.    10.5      35.5
## 5 Dry                  After Experiment       371.     9.21     26.7
## 6 Dry                  After Rehydration      371.     9.23     22.7

SMI

ggplot() + 
  geom_point(data = rehydrat_dat,
             aes(x = day,
                 y = SMI, 
                 color = humidity_tmt_percent
                 ), 
             size = 1, 
             alpha = 0.6) + 
  geom_line(data = rehydrat_dat,
            aes(x = day,
                y = SMI,
                group = individual_ID,
                color = humidity_tmt_percent),
            alpha = 0.2) +
  geom_line(data = hydrat_means, # new data
            aes(x = day,
                y = SMI_mean,
                group = humidity_tmt_percent,
                color = humidity_tmt_percent),
            size = 1.6,
            alpha = 1) +
  theme_classic() + 
  scale_color_brewer(palette = "Set2",
                     name = "Humidity Treatment") +
  scale_x_discrete(labels = c("Before\nExperiment",
                              "After\nExperiment",
                              "After\nRehydration")) +
  xlab("") + 
  xlab("") + 
  ylab("Body Condition (g)") + 
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 18),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 14),
        legend.text.align = 0,
        legend.position = "bottom"
        ) -> rehab_SMI_fig
rehab_SMI_fig

# export figure
#ggsave(filename = "rehab_SMI_fig.jpeg",
 #      plot = rehab_SMI_fig,
  #     path = "./final_figures",
   #    device = "jpeg",
    #   dpi = 1200,
     #  width = 5, height = 4)

Osmolality

first, make a list of all the IDs that have a post-rehab osmolality measurement, since this has a lot of missing data

rehab_osmols <- rehydrat_dat %>%
  dplyr::filter(day == "After Rehydration") %>%
  dplyr::filter(complete.cases(osmolality_mmol_kg))
rehydrat_dat %>% 
  dplyr::filter(individual_ID %in% rehab_osmols$individual_ID) %>%
  ggplot(data = .) + 
  geom_point(aes(x = day,
                 y = osmolality_mmol_kg, 
                 color = humidity_tmt_percent
                 ), 
             size = 1, 
             alpha = 0.6) +
  geom_line(aes(x = day,
                y = osmolality_mmol_kg,
                group = individual_ID,
                color = humidity_tmt_percent),
            alpha = 0.2) +
  geom_line(data = hydrat_means, # new data
            aes(x = day,
                y = osml_mean,
                group = humidity_tmt_percent,
                color = humidity_tmt_percent),
            size = 1.6,
            alpha = 1) +
  theme_classic() + 
  scale_color_brewer(palette = "Set2",
                     name = "Humidity Treatment") +
  scale_x_discrete(labels = c("Before\nExperiment",
                              "After\nExperiment",
                              "After\nRehydration")) +
  xlab("") + 
  xlab("") + 
  ylab("Osmolality (mmol/kg)") + 
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 18),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 14),
        legend.text.align = 0,
        legend.position = "bottom"
  ) -> rehab_osml_fig
rehab_osml_fig

# export figure
#ggsave(filename = "rehab_osml_fig.jpeg",
 #      plot = rehab_osml_fig,
  #     path = "./final_figures",
   #    device = "jpeg",
    #   dpi = 1200,
     #  width = 5, height = 4)

Hematocrit

first, make a list of all the IDs that have all three measurements:

rehab_hct <- rehydrat_dat %>%
  dplyr::filter(day == "After Rehydration") %>%
  dplyr::filter(complete.cases(hematocrit_percent))
rehydrat_dat %>% 
  dplyr::filter(individual_ID %in% rehab_osmols$individual_ID) %>%
  ggplot(data = .) + 
  geom_point(aes(x = day,
                 y = hematocrit_percent, 
                 color = humidity_tmt_percent
                 ), 
             size = 1, 
             alpha = 0.6) +
  geom_line(aes(x = day,
                y = hematocrit_percent,
                group = individual_ID,
                color = humidity_tmt_percent),
            alpha = 0.2) +
  geom_line(data = hydrat_means, # new data
            aes(x = day,
                y = hct_mean,
                group = humidity_tmt_percent,
                color = humidity_tmt_percent),
            size = 1.6,
            alpha = 1) +
  theme_classic() + 
  scale_color_brewer(palette = "Set2",
                     name = "Humidity Treatment") +
  scale_x_discrete(labels = c("Before\nExperiment",
                              "After\nExperiment",
                              "After\nRehydration")) +
  xlab("") + 
  ylab("Hematocrit (%)") + 
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 18),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 14),
        legend.text.align = 0,
        legend.position = "bottom"
  ) -> rehab_hct_fig
rehab_hct_fig

# export figure
#ggsave(filename = "rehab_hct_fig.jpeg",
 #      plot = rehab_hct_fig,
  #     path = "./final_figures",
   #    device = "jpeg",
    #   dpi = 1200,
     #  width = 5, height = 4)

Multi-Figure

ggarrange(rehab_osml_fig, rehab_SMI_fig, rehab_hct_fig,
          ncol = 1, nrow = 3,
          common.legend = TRUE,
          legend = "bottom"
          ) -> rehab_multi_fig
rehab_multi_fig

# export figure
#ggsave(filename = "rehab_multi_fig.jpeg",
 #      plot = rehab_multi_fig,
  #     path = "./final_figures",
   #    device = "jpeg",
    #   dpi = 1200,
     #  width = 6, height = 12)